Large scale structure simulations of inhomogeneous LTB void models 
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We perform numerical simulations of large scale structure evolution in an inhomogeneous 
Lemaitre-Tolman-Bondi (LTB) model of the Universe. We follow the gravitational collapse of a 
large underdense region (a void) in an otherwise fiat matter-dominated Einstein-deSitter model. 
We observe how the (background) density contrast at the centre of the void grows to be of order 
one, and show that the density and velocity profiles follow the exact non-linear LTB solution to the 
full Einstein equations for all but the most extreme voids. This result seems to contradict previous 
claims that fully relativistic codes are needed to properly handle the non-linear evolution of large 
scale structures, and that local Newtonian dynamics with an explicit expansion term is not ade- 
quate. We also find that the (local) matter density contrast grows with the scale factor in a way 
analogous to that of an open universe with a value of the matter density Qmir) corresponding to 
the appropriate location within the void. 



PACS numbers: 98.80.Cq 

I. INTRODUCTION 

Distant supernovae appear dimmer than expected in 
a purely matter-dominated homogeneous and isotropic 
FRW universe. The currently favoured explanation of 
this dimming is the late time acceleration of the universe 
due to an energy component that acts like a repulsive 
force. The nature of the so-called Dark Energy responsi- 
ble for the apparent acceleration is completely unknown. 
Observations seem to suggest that it is similar to Ein- 
stein's cosmological constant, but there is inconclusive 
evidence pQ. In the meantime, our realization that the 
universe around us is far from homogeneous, since there 
are large superclusters and huge voids across our largest 
galaxy catalogs [2] has triggered the study of alterna- 
tives to this mysterious energy. Since the end of the 
nineties it has been suggested by various groups [31 3] 
that an isotropic but inhomogeneous Lemaitre-Tolman- 
Bondi universe could also induce an apparent dimming 
of the light of distant supernovae, in this case due to lo- 
cal spatial gradients in the expansion rate and matter 
density, rather than due to late time acceleration. There 
is nothing wrong or inconsistent, apart from philosophi- 
cal prejudices, with the possibility that we live close to 
the centre of a gigaparsec-sized void. Such a supervoid 
may indeed have been observed as the CMB cold spot [5] 
and somewhat smaller voids have been seen in the local 
galaxy distribution [5J [7] . If a local void had the size and 
depth of a void responsible for the cold spot, i.e. r ~ 2 
Gpc and Ojvf ~ 0.2 within a flat Einstein-de Sitter uni- 
verse, it would be consistent with local observations [51(5], 
and could account for the supernovae dimming, together 
with the observed baryon acoustic oscillations and CMB 
acoustic peaks, the age of the universe, local rate of ex- 
pansion, etc. [il [rUrtT3] . 

In order to make contact with large scale structure 
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observations of the matter distribution, large numerical 
simulations are usually performed, where a very specific 
initial condition is assumed for the primordial spectrum 
of inhomogeneities and the evolution is done solving the 
Newtonian dynamics numerically. In most cases, the 
matter content is just cold dark matter falling into gravi- 
tational potential wells set in by inflation, although some 
simulations have included also baryons as well as hot dark 
mater, neutrinos, radiation, and astrophysical feedbacks. 

This conceptually simple recipe yields results which 
are in good agreement with the matter distribution we 
observe in the sky on large scales, and can be used to 
constrain our model of the Universe and determine some 
of the parameters of the Standard Model of Cosmology. 
However, some have argued (see e. g. [TS] and references 
therein) that the late stages of gravitational collapse, and 
structure formation in the universe require a fully rela- 
tivistic numerical description in order to capture specific 
signatures of the strong non-linear dynamics of general 
relativity, and make a correct treatment of features even 
with sizes comparable to the Hubble radius, see also [T71 - 

EH]. 

In this paper we have tested the validity of the New- 
tonian approximation for structure formation in the con- 
text of an inhomogeneous model whose fully non-linear 
dynamics can be solved exactly using the Einstein equa- 
tions [3, 4J. We start with an Einstein-deSitter model at 
a high redshift, and we test it under various initial condi- 
tions: i) We first include just a large void of fixed size and 
a small initial amplitude, and we follow the non-linear 
growth of the void's depth and size; ii) We then add Cold 
Dark Matter (CDM) with a Gaussian random field dis- 
tribution based on inflation (n s = 1 and as = 0.9), with 
the weighted matter-baryon transfer function included 
to account for the Baryon Acoustic Oscillations (with 
/g as = 0.14 and Q to t = 1 consistent with WMAP-7yr), 
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and follow the growth of both the void and the matter 
power spectrum. We have confirmed that our numerical 
simulations follow the exact solution of the LTB back- 
ground Einstein equations at all scales (radii), except for 
very extreme cases, where shell crossing occur in mod- 
els with large scale structure fluctuations (shell crossing 
does not occur though if we take a pure void) . This seem 
to suggest that the Newtonian approximation for grav- 
itational collapse is perfectly valid, even for gigaparsec- 
sized voids as empty as J1m = 0.02 at the centre at z = 0, 
which corresponds to density contrasts of order 1 with re- 
spect to the asymptotic EdS model. We obtain very good 
matches to both the density and the velocity profile for 
matter moving in such LTB backgrounds. We also check 
that the non-linear evolution gives rise to well differenti- 
ated Hubble rates along the line of sight and transverse 
directions, in perfect agreement with the exact relativis- 
tic solutions. Moreover, we can follow the evolution of 
the matter density contrast as a function of the scale fac- 
tor and find that it evolves as one would expect for an 
open universe with the value of f2 m (t) corresponding to 
the local position within the void. 



II. LEMAITRE-TOLMAN-BONDI VOID 
MODELS 

The Lemaitre-Tolman-Bondi model describes general 
spherically symmetric space-times and can be used as a 
toy model for describing large voids in the universe. The 
metric is given by 



ds 2 = -dt 2 



A' 2 (r,t)dr 2 
1 - k(r) 



A 2 (r,t) (d9 2 



with a spherically symmetric matter source with negligi- 
ble pressure, T% = — pm{t, t) Sq 5°. Since we have dif- 
ferent radial and angular scale factors, we also define a 
transverse and longitudinal Hubble rates as Ht = A/A, 
and Hl = A' /A' , where dots and primes denote dt and 
d r , respectively. Integrating the Einstein equations for 
this metric one finds the r-dependent transverse Hubble 
rate 



( Mr) 
\A(r,t) 



where we have fixed the gauge by setting Ao(r) = r and 
Qi((r) = 1 — fljif(r). For fixed r the above equation is 
equivalent to the Friedmann equation, and has an exact 
parametric solution, see Ref. [I]. 

In general, LTB models are uniquely specified by the 
two functions Hq(t) and f2jy-(r), but to test them against 
data we have to parameterize the functions, to reduce the 
degrees of freedom to a discrete set of parameters. For 
simplicity in this paper we will use the constrained GBH 
model [4] to describe the void profile. First of all, it uses 
a minimum set of parameters to make a simple void pro- 
file, and secondly, we impose that the time to Big Bang 



should be constant. We have made this choice, because 
models with an inhomogenous Big Bang would contain 
a mixture of growing and decaying modes, and conse- 
quently the void would not disappear at high redshift, 
making them incompatible with the Standard Big Bang 
scenario [20] ■ If we only consider constrained LTB mod- 
els, then at high redshifts and/or at large distances the 
central void is reduced to an insignificant perturbation 
in an otherwise homogeneous universe described by an 
FRW metric, and physical results for the early universe 
derived for FRW space-times still hold, even though we 
are considering an LTB space-time. The second condi- 
tion gives a relation between Ho(r) and fijvf (f), and hence 
constrain the models to one free function, and a propor- 
tionality constant describing the overall expansion rate. 
Our chosen model is thus given by [U [TP] 
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where we have assumed that space is asymptotically flat, 
fi(oo) = 1. The model has then only four free parame- 
ters: The overall expansion rate H , the underdensity at 
the centre of the void fli n , the size of the void ro, and the 
transition width of the void profile Ar. For more details 
on the model see Ref. [4]. 



III. LINEAR PERTURBATION THEORY 

We still do not have a complete linear perturbation the- 
ory for LTB models. The main difficulty is that since the 
background is inhomogeneous we cannot split the per- 
turbations into independent equations for the scalar, vec- 
tor and tensor modes. In LTB models the equations for 
these modes appear as coupled partial differential equa- 
tions [201 IS]- In particular, the scalar modes couple to 
the tensor shear modes at first order, which act as source 
for the scalar mode via the background shear. However, 
in the case that the latter is small, like in the models 
we have been describing in our previous works we 
can ignore this source and solve exactly the perturbation 
equation for the scalar mode $, which in the absence 
of anisotropic matter stresses is equal to the curvature 
mode ^ . In this approximation the equation becomes 



<£(r,i) + 4H T (r,t)$(r,t) 



with the exact solution 



2k(r) 
A 2 (r,t) 



*(r,*) = J (1) 



$(r, t) = $ (r,0)aFi 



1,2, 



A{r,t) 



(2) 

We note that, strictly speaking, this solution is only ex- 
act when ignoring the tensor coupling, and considering 
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angular transverse modes, but turns out to be a very 
good approximation. In that same approximation (neg- 
ligible background shear) , the density contrast of matter 
is proportional to the scalar metric perturbation, 



S(r,t) = S (r) 



A(r,t) $(r,t) 
r $ (r,0)' 



(3) 



where Sq(t) 7 up to a normalization factor, can be deter- 
mined under the assumption that the small scale mat- 
ter perturbations in the early universe decouple from the 
void, giving 5o(r) oc r / A{r,t car i y ). It is this function 
which we will try to compare with the simulations de- 
scribed in the next section. 



Name 


2start 




Ar/r 


^particles 


Comments 


H 


24 


0.25 


0.3 


960 3 


High res sim 


V 


49 


0.25 


0.3 


512 3 


Void alone 


524 


24 


0.25 


0.3 


512 3 


Void + matter 


549 


49 


0.25 


0.3 


512 3 


Void + matter 


599 


99 


0.25 


0.3 


512 3 


Void + matter 


5Q125 


49 


0.125 


0.3 


512 3 


Void + matter 


5^063 


49 


0.0625 


0.3 


512 3 


Void + matter 


5^021 


199 


0.0208 


0.3 


512 3 


Void + matter 


5A01 


49 


0.125 


0.1 


512 3 


Void + matter 


5A05 


49 


0.125 


0.5 


512 3 


Void + matter 


C 


49 


0.25 


0.3 


768 3 


L=3600 Mpc h' 1 



IV. NUMERICAL SIMULATIONS 

To test the validity of N-Body codes in describing gi- 
gaparsec sized voids, and to follow the evolution and for- 
mation of structure in such models, we have modified the 
2LPT initial condition generator [22] to set up an N-Body 
simulation of a void for the Gadget2 code [53] where the 
displacements and velocities of the particles are found 
using second order Lagrangian perturbation theory |22j . 
Starting with a standard transfer function for the total 
matter content in a flat Einstein-deSitter model we con- 
struct initial conditions for the gravitational potential in 
fc-space . Then we find the gravitational potential of a 
void using the analytical solution, by interpolating the 
density out on the particle grid, and then Fourier trans- 
forming it. Now that the total potential $k = $ k + 
is known, the 2LPT code proceeds unchanged from the 
original version. Once the initial conditions have been 
set up we use the public domain version of the Gadget2 
code in pure tree-mode to run the simulation (see table 
| for an overview of the simulations) [28] . 

It is not evident that N-body simulations can be used 
to describe large scale LTB models, and therefore a signif- 
icant effort has gone into validating that indeed we repro- 
duce the expected theoretical behaviour. To test the code 
we have used different starting redshifts (z start =24, 49, 
99, and 199) to check explicitly that the code is started 
at high enough redshifts, such that the displacements of 
the particles are much smaller than the inter-particle dis- 
tance, and that the void can be treated as a linear per- 
turbation, which at first order does not interact with the 
small scale fluctuations from the power spectrum. We 
have used different resolutions (simulations 524 and H - 
see table |l| to test that the cosmological large scale struc- 
ture is adequately resolved, we have tested that the void 
does not interact too much with mirror images of itself by 
changing the physical box size from L=2400 to £=3600 
Mpc h^ 1 (simulations 549 and £), and we have checked 
that to first order the small scale fluctuations do not 
back react significantly on the void, by running with and 
without matter perturbations (simulations 549 and V). 
Apart from the numerical tests, we have simulated a rep- 
resentative set of realistic void models varying the tran- 



TABLE I: Overview of the simulations. All have been per- 
formed with a void of radius ro = 1100 Mpc = 473 Mpc ft -1 , 
and with an asymptotic Hubble parameter of hoo = 0.43. The 
standard box size is L=2400 Mpc h~ , and the particle mass 
is Mpart= 2.8 x IO^Mq/i" 1 (M part = 4.3 x lO^Af©^ 1 for 
H). Everywhere we have used a smoothing length of 56 kpc 
ft -1 (except for H, where it has been appropriately rescaled). 

sition length Ar/ro and central underdensity fl ln (see ta- 
ble [I]). The majority of the simulations use a GBH model 
with Sl in = 0.25, and Ar/r = 0.3, but we have also 
run other simulations with f2j n = 0.125, f2i n = 0.0625, 
il in = 0.0208, and Ar/r = 0.1 and Ar/r Q = 0.5. 



V. ANALYSIS AND RESULTS 

The results in this paper show the concordance be- 
tween the simulations and the theoretical predictions. In 
order to check this we use the highest resolution simula- 
tion H. as our reference model and the other simulations 
to test the limits of the validity of N-Body simulations for 
describing LTB models. We have subjected our simula- 
tion to three different tests, confronting the density pro- 
file, the Hubble parameter profile (Ht and Hl) and the 
density contrast evolution with the corresponding theo- 
retical predictions. 

A. Distances and redshifts in LTB models. 

Gadget2 is designed to perform simulations of FRW 
universes, and one needs to associate the comoving ra- 
dial coordinates in both models. Since we have analyzed 
the data from Gadget2 snapshots, that is, positions and 
velocities of particles are "measured" at constant cosmic 
time, and all our observables are quantities calculated 
in thin spherical shells, this identification must be done 
through the proper radial distance, calculated in both 
cases as 

d p (r,t)=a(t)r FRW = \' dr. (4) 

Jo V 1 _ k V) 
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FIG. 2: Comparison of the density profile of the H simulation 
at different redshifts with the theoretical curves, as a function 
of comoving distance tfrw = (1 + z)d p in Mpc. 



Proper distance [Mpc h -1 ] 



FIG. 1: The projected matter distribution at z — averaged 
over a 175 Mpc slice centered on the void of the 960 3 simu- 
lation. Notice how near the centre of the void, not only the 
density is lower, but also there is significantly less structure 
than outside the void. The characteristic void size ro = 473 
Mpc is indicated by the thin circle. 



When the curvature factor (1 — fc(r)) -1 / 2 is roughly 1, 
which is the case in the models under study, one can 
approximate 

d p (r, t) = a(t) r F Rw - ^(jltb, *) (5) 

for most redshifts. Similarly, when interpreting the re- 
sults, it is important to remember that while the proper 
cosmological time in the two metrics can readily be iden- 
tified, the redshift at equal times are different, i.e. for 
^frw = £ltb the zfrw an d zltb are different. It is 
important to emphasize that since we are considering a 
constrained-GBH LTB model, the time to Big Bang is 
homogenous and thus all times at each radial distance 
are the same, so each particle in the simulation has a 
time given by the code: £frw = *ltb- 

B. Density profile 

We first compare the theoretical density profile of a 
GBH universe having the desired parameters with the 
corresponding profile obtained from the simulation. The 
density field is calculated by interpolating each particle 
in the box to a grid using a 2 nd order triangular-shaped- 
cloud (TSC) technique [21] (see fig. [TJ; then the simula- 
tion box is divided into different spherical bins, and we 
calculate the average density in each of them thus ob- 
taining the density as a function of the proper distance 
dp (see eq. [41) . Due to the presence of non-linear inhomo- 
geneities, the error in the determination of the density 



profile cannot be directly obtained as the r.m.s. in each 
bin, and the error bars displayed in the figures have been 
calculated as the r.m.s. in the analogous V simulation 
without CDM perturbations. The reference simulation H 
shows an excellent agreement between theory and simula- 
tion (see fig. [2]) , except near the centre of the void, where 
the particle distribution is undersampled and shot noise 
dominated. 

In fig. [3] we show the density profile for an extended 
set of models. For most models the simulations are in 
excellent concordance with the theory, though for two 
extremal cases, namely the emptiest void 517021, and 
the void with the steepest transition SOA01 we find sig- 
nificant deviations. For SS7A01 the discrepancy is not 
severe, and only present in the density profile. We spec- 
ulate that this could be due to under resolution of the 
transition length or possibly due to the small-scale per- 
turbations interacting with the large scale void, given 
that the transition length is only Ar = 47.3 Mpc h -1 . 



C. Rates of expansion 

The radial velocity profile can be used to compare 
against the theoretical predictions for and H^. The 
rate of change in the proper distance d p /d p computed in 
the rest-frame of the matter should match each other in 
the FRW and LTB metric, if the simulations are a valid 
description of the LTB model. In the LTB metric matter 
is at rest, and keeps the same comoving coordinate, while 
in the FRW metric there are systematic radial motion, 
and We have that 

^dp RW = — [ar matter ] = d p [{v r )/r + H^] , (6) 

which can be directly compared to the theoretical LTB 
result calculated taking the derivative of the r.h.s. of 
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FIG. 3: Density profiles for different values of Sli n (top panel) 
and Ar/ro (lower panel) in comparison with the correspond- 
ing theoretical profiles. All curves are plotted at redshift 
z = 0. 



FIG. 4: The velocity profile for the simulation H (top panel) 
and the Ht and Hl profiles for different values of fli n (lower 
panel). In both cases, the theoretical profiles are shown with 
dotted lines. All curves are plotted at redshift z = 0. 



eq. Q. (v r ) is calculated as the average radial velocity v r 
of the particles sampled in spherical bins. In the upper 
panel of fig. [4] we see how the theoretical radial velocities 
(calculated from a _1 [d p — -ffoodp] ) match the data from 
H. We have found that d p /d p aproximate Ht very well 
(see eq. [5]), possibly because the denominator in eq. Q, 
(1 — ^(rjj 1 / 2 , being time independent, cancel in the ra- 
tio dp I dp. Using this approximation in the lower panel of 
fig. [4] we compare a range of models to theory. Again, the 
difference with the theoretical graph found near d p = 
is understandable, we are shot noise dominated, and fur- 
thermore the matter perturbations displace the centre of 
the void slightly, while at the same time we have a formal 
singularity at r — when calculating (v r )/r. From Ht 
we can extract Hl straightforwardly as : 

Hl = §=Ht + ^H t , (7) 

which is just Hl — Ht + r H' T at z = 0, using Aq(v) = r. 
We stress though, that at z = this is a derived param- 
eter, and not independent of We find that all 
but the emptiest model 6>f2021 match well the theoret- 



ical predictions. For <Sf2021 the velocity is consistently 
higher (and the density lower) inside the void compared 
to theoretical predictions, and a density spike is building 
up near the edge of the void. This could be due to the 
very lower density, but it may also be a consequence of 
the very high starting redshift (z B tart = 199), that was 
necessary to keep the perturbations linear and the parti- 
cle displacements acceptable in the initial condition. 



D. Density contrast 

Another interesting observable to study is the evolu- 
tion of the density contrast as a function of redshift, 
8(z) = ((fi(z) — p)/p). Being (random) fluctuations, we 
calculate it by finding the r.m.s. of S(z) in spherical bins 
as a function of proper distance. The errors in the deter- 
mination of S were calculated as the standard deviation of 
the values of 6 calculated in the 8 octants of each spher- 
ical bin. The results for the simulation 549 can be seen 
in fig. [5j where we compare the density contrast, calcu- 
lated at a fixed comoving distance rp^w = (1 + z) d p , 
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FIG. 5: Density contrast evolution inside the void at comov- 
ing distance tfrw = + d v = 280 Mpc, for 549, in compar- 
ison with the theoretical prediction from perturbation theory 
(full line). We also compare the LTB growth of density per- 
turbations with that of Open CDM (dotted line) and ACDM 
(dashed line) . The theoretical curves were normalized to have 
the same slope asymptotically in the past, as apRw — > 0. Note 
that, even though the horizontal axis reads 1/(1 + ,zfrw), this 
zfrw only determines the cosmic time t, since the density 
contrast was calculated at a fixed comoving distance, and not 
in the lightcone. 



as a function of time (expressed in terms of rcdshift), 
with the predicted one within the simplified linear per- 
turbation theory in LTB described by eq. ([3]). We also 
include, for comparison, the density contrast growth for 
an open universe, with Qm = 0.25, and a ACDM, with 
£Im = 0.25 and fl\ = 0.75. It is interesting to note 
that the data agree well, within error bars, with both 
the theoretical prediction in the LTB model and in the 
concordance ACDM, while they differ significantly from 
an open universe with the same matter density. 

In order to better understand these differences, we have 
also studied the evolution of the density contrast at sev- 
eral distances from the center of the void. The results are 
shown in fig. [6] Two clearly different zones can be dis- 
tinguished: while the growth is proportional to (1 + z) _1 
for large comoving distance, i.e. 51 ~ 1 outside the void 
and apRw = (1 + z ) _1 'j the growth is significantly slower 
(as would occur in an open FRW universe), for small 
distances. 



VI. CONCLUSIONS 

We have studied for the first time the non-linear evo- 
lution of structure formation in large-void LTB models 
within an asymptotic Einstein-deSitter (EdS) universe. 
By initiating large N-body simulations at high redshifts, 
we have been able to follow the non-linear gravitational 
collapse of matter structures in the presence of an un- 
derdense void that starts with a density contrast of or- 



0.6 




1/(1 

FIG. 6: Density contrast evolution at different fixed comov- 
ing distances for 549 simulation, as a function of the FRW 
scale factor, like in Fig. 5. It is easy to distinguish between 
the contrast growth in a background with Q, ~ I at large 
distances and with a lower Q near the void center. 



der S m ~ 10~ 3 at photon decoupling (where the mat- 
ter perturbations have S m ~ 10~ 5 ). We find that us- 
ing a standard N-body code, the nonlinear growth of the 
void underdensity follows the exact analytical solution 
of the Einstein equations, even for very deep voids with 
CIm = 0.06 at the center, and thus with density contrasts 
of order one with respect to the asymptotic EdS universe. 
Moreover, the transverse and longitudinal rates of expan- 
sion agree with the theoretical expectations, giving us 
confidence that the simulations are tracing the full non- 
linear gravitational collapse in this non-perturbative LTB 
background. This is furthermore evidence that N-body 
codes give a credible and precise description of the stan- 
dard ACDM model, where the voids are of much lesser 
size, and no general relativistic corrections are needed to 
describe the large scale evolution. 

We have also studied the evolution of the matter den- 
sity contrast in such a non-trivial background, and found 
an analytical solution to the approximate equations for 
the growth of perturbations in the limit of negligible 
background shear, and shown that the numerical and 
analytical results are in good agreement. Moreover, the 
comparison with OCDM and ACDM shows that the den- 
sity contrast growth for our LTB models is very close 
within errors to that of the concordance ACDM models 
suggested by WMAP-7yr Q]. 

From our non-linear LTB N-body simulations we can 
extract predictions for observations of large scale struc- 
ture, via the two-point angular correlation function, the 
angular power spectra, the growth of structure, and the 
density contrast. Our models therefore give the possi- 
bility of using both current and future observations of 
the large scale structure (such as DES [55], EUCLID |H] 
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and PAU [27 ), to constrain LTB models that already 
provides a viable fit to current observations of the geom- 
etry of the Universe. 
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